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Abstract 

In this paper we introduce a class of stochastic population models based on "patch dynamics." 
The size of the patch may be varied, and this allows one to quantify the departures of these 
stochastic models from various mean field theories, which are generally valid as the patch size 
becomes very large. These models may be used to formulate a broad range of biological processes 
in both spatial and non-spatial contexts. Here, we concentrate on two-species competition. We 
present both a mathematical analysis of the patch model, in which we derive the precise form of 
the competition mean field equations (and their first order corrections in the non-spatial case), and 
simulation results. These mean field equations differ, in some important ways, from those which 
are normally written down on phenomenological grounds. Our general conclusion is that mean 
field theory is more robust for spatial models than for a single isolated patch. This is due to the 
dilution of stochastic effects in a spatial setting resulting from repeated rescue events mediated by 
inter-patch diffusion. However, discrete effects due to modest patch sizes lead to striking deviations 
from mean field theory even in a spatial setting. 

PACS numbers: 87.23.Cc,02.50.Ey,05.40.-a 
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I. INTRODUCTION 



Traditional theoretical ecology, in which the time evolution of population densities is 
described by differential equations, has a long history For a single species the 

simplest form of the governing equation is assumed to take the form dN/dt = ^{N)N, 
where $(A^) describes the growth of the population. A common choice when modeling this 
growth is to take ^{N) = r(l — N/K), where r and K are two constants. By analogy, when 
describing the interaction of two species, it is natural to postulate that the populations 
of the two species, Ni and A''2, change according to dNi/dt = f{Ni,N2) and dN2/dt = 
g{Ni,N2). The functions / and g are chosen according to whether the interactions are 
purely competitive, predator-prey like or include other effects such as co-operation. We 
will refer to descriptions of this kind as population-based; they are arrived at without the 
need for a detailed knowledge of the interaction between individuals and rely instead on 
assuming that the terms which arise in the governing equations represent the net effects of 
individual interactions in some generic way. Equations of this kind play such a central role 
in population biology, that many subsequent elaborations of the theory have taken them as 
the starting point: spatial variation is introduced by adding a drift term V'^N^ (« = 1, 2) 
to the right-hand side of the ath equation, and the models are sometimes interpreted as 
referring to individuals by assuming that the functions / and g also describe interactions at 
the level of the individual. 

In the last decade or so, an alternative approach to that of classical theoretical ecology 
described above has been developed. This involves abandoning the traditional population- 
level description in favor of an individual-based description in which explicit rules governing 
the interaction of individuals with each other and with the environment are given. The 
popularity of these individual-based models (IBMs) is undoubtedly due to the continuing 
increase in the availability of powerful computers, but they also have other attractive fea- 
tures, such as the ability to directly model individual attributes. At this point we should 
stress we are assuming that the individuals of a given species in our model are identical, 
and thus the term IBM should not be confused with agent based models which are often 
designed to study the ecological effects of behavioral and physiological variation among in- 
dividuals. A better term might be ILM (individual level model), but the term IBM has wide 
usage, and so we will use it here. In this paper we will be concerned with theoretical issues 
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which relate to the connection between models defined at the individual level and those at 
the population level. Thus, the individuals in our models will be identical within a given 
species. The relation between population-level and individual-level descriptions has been 
a focus of discussion within the theoretical ecology community for some time jj, 0]. Some 
regard the nature of the population-based models as obvious and either write them down 
without comment, or derive mean-field equations by making an assumption of homogeneous 
mixing of the populations 0, Q] . However, there is also some recognition that the situation 
may be more complicated than this 0, 0] , and that the transition to a partial differential 
equation required for a spatial description, from the ordinary differential equation obtained 
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by using mean-field theory, may not be as simple as just adding the term V^A^^i 

From a statistical physics perspective it is natural to expect that fiuctuations may play 
an important role in these systems and lead to important differences between microscopic 
(IBM) and macroscopic (population-level) descriptions. A formal analysis of such issues has 
been presented for simple birth/death processes and annihilation reactions Q| using the 
language of field theory and the renormalization group. Here, we take a different approach, 
using van Kampen's system size expansion, in order to probe the connection between more 
complex IBM's and mean field theories. This has some advantages, for instance it is not 
necessary using this technique to first construct a corresponding field theory, and then 
extract the mean-field equations and the Gaussian fiuctuations about them to model the 
macroscopic behavior. There is also less of an emphasis of the role of phase transitions, 
which are not usually of prime interest in ecological models. 

We would like to stress from the outset that our use of "mean field theory" or "mean 
field model" is in the spirit of statistical mechanics. By "mean field" we mean the neglect 
of correlations between degrees of freedom, allowing one to write the mean of the product 
of two stochastic variables as the product of their means. Thus, a spatial model in which 
statistical correlations are neglected is still a "mean field model" within this usage. This is in 
contrast to some papers in population biology in which "mean field" is reserved exclusively 
to describe non-spatial models. 

In this paper, the models which we study will be defined only at the level of direct 
interaction among individuals and in terms of local properties such as birth, death and 
migration rates. The population-based properties of the model will then be derived within 
a well-defined approximation scheme. We will show that, beginning with reasonable models 
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at the individual level, the corresponding population-level models are similar in structure 
to those that we would naively write down, but have important differences. For instance, 
small scale diffusion may not simply translate into V'^N^ terms and the parameters defined 
at the individual-based level may not map directly into their equivalents at the population 
level. 

The individual interactions may be naturally introduced using a "patch model." From 
a biological perspective, a patch can be thought of as a small spatial region within which 
interactions between individuals occur. The patch is assumed to be sufficiently small that 
there are no spatial effects. In other words, there is complete mixing, and all individuals 
have the same chance of potentially interacting with each other. In the non-spatial version 
of the model this simply amounts to stating that the probability of any individual dying in 
unit time should be proportional to the density of individuals existing at that time, and for 
processes which involve two individuals, the probability involved should be that found when 
drawing two of these types of individuals at random from a patch which contains all the 
individuals in the system. Not all constituents of the patch will correspond to individuals; 
some will correspond to empty sites in the spatial version of the model. A more detailed 
specification is given in Sectional Patch models such as this have been used in many areas 
of science Q], and often go under the name of "urn models." Two early examples were the 
Ehrenfest urn, which was used to discuss the foundations of statistical mechanics, and the 
P61ya un, wKicK was on„„a,.v devised to de.e.be ooMagic. disease. Q. both ease, 
the problem of interest can be mapped onto an urn which contains balls of different colors, 
say black and white, which are drawn randomly one at a time. In the case of the Ehrenfest 
urn, each time a ball is drawn it is replaced by one of the other color. For the Polya urn, 
the chosen ball is replaced together with an extra ball of the same color. The relation to 
the modeling of contagious diseases should be clear: each occurrence of a particular color 
increases the probability of further occurrences. Since the introduction of these particular 
urn models, the notion has been generalized considerably. For example, both models fall 
into the class of urn models which have the properties that if a white ball is drawn, it 
is replaced together with a white balls and b black ones, and if a black ball is drawn, it is 
replaced together with c white balls and d black ones. A further generalization is to consider 
r different colors, with the drawn ball of color i being replaced together with aij balls of 
color j = 1, . . . ,r) It is also clear that two, or more, balls can be drawn together 
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or, as we will have occasion to assume in this paper, two balls could be drawn for a fraction 
fi of the time and one ball for a fraction (1 — /i) of the time. 

Urn models are concrete realizations of stochastic processes with probabilities which 
depend only on the instantaneous state of the system, that is, Markov processes. They 
have proved useful in several areas of the biological sciences. Perhaps the most obvious 
application is in population genetics; implicitly in the early work of Fisher and Wright , 
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and explicitly in later developments [19|, |20|, |21|, |22, |23|, |2J, |25[. However they have also 
)een used in a number of other areas such as the study of radioactive particles in animals 



2^! 21> the study of patterns of vegetation [2^, models of interaction between species 
]) 31, 32 1 and metapopulation models 3^. The last three applications are closest to the 
ones discussed in this paper, but differ in that the balls represent forest or grassland in the 
first case, species in the second case and colonies in the third case, rather than individuals 
as in the present paper. 

Since urn, or patch, models are representations of Markov processes, the continuous 
time version of their dynamics may be described using a master equation. The use of 



master equations is familiar to physicists, but they are sti 
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biological community (but see above references and Refs. 
process has been formulated in this way, we may use standard techniques to take the mean 
field limit and so obtain the corresponding population-based equation. Much of this paper 
will be taken up with a comparison of results obtained from the full model (averaged over 
many realizations) and the mean-field results. Our approach will be particularly useful in 
distinguishing situations where mean-field theory is a good approximation to the full theory 
and situations where it is not. We will mainly concentrate on competitive interactions in 
systems with one or two species, but our method applies equally well to predator-prey or 
epidemic models and multispecies communities. 

The outline of the paper is as follows. We first motivate and define the rules for the non- 
spatial patch models in Section m starting with a single species model, and then generalizing 
to a two-species competition model. The full stochastic nature of the models is described 
using master equations, and the associated mean field models are derived. In the following 
section we go one step beyond mean field theory and derive the dynamics of the Gaussian 
fluctuations about the mean field solutions. This is equivalent to not only following the mean 
position of the probability distributions over time, but also describing their broadening. 
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In Section IIVI we present numerical simulations of the fully stochastic non-spatial models, 
for both one and two species. We compare our results to numerical integration of both 
the mean field equations, and the improved models with Gaussian fluctuations included. 
We find that the agreement between mean field theory and simulations is excellent for 
larger patches, as is to be expected, but that the breakdown of mean field theory occurs 
precipitously below critical patch sizes which are still fairly large. We proceed in Section IVl 
to generalize our patch model formulation to spatially explicit population dynamics of two 
species competition. From the master equation we derive the spatial mean field equations. 
We note important differences between these equations and "intuitive" versions which have 
appeared in the literature. These mean field equations are tested against simulations of 
the fully stochastic models in Section IVll In this introductory paper we are unable to 
give a comprehensive analysis of the two species model. Instead we present two interesting 
scenarios, and note the successes and failures of the mean field equations, which are strongly 
dependent on patch size and the presence/absence of interspecific competition. We end the 
paper with our conclusions along with a discussion of future directions. Two appendices 
contain technical details. The first concerns the system size expansion for two species and 
the second the formalism for spatial systems. 



II. BASIC FORMALISM AND THE NON-SPATIAL MODEL 

In this section we will introduce the essential features of our approach by formulating an 
individual based stochastic model of competition between two species. We will then show 
that in the limit of large population sizes, the time evolution of this model reduces to the 
well known differential equations describing the population growth of two competing species. 

The two species will be labeled A and B. To motivate the approach we will adopt, suppose 
that we model the interactions of A and B individuals in an area of land by subdividing 
it into plots of equal area. The plot sizes are chosen so that each one either contains 
one A individual, or one B individual, or neither an A nor a B. We will call the latter an 
empty site and label it by E. In a spatial version of the model we would give rules for how 
A, B and E interact, specify birth and death rates for A and B, and allow them to move to 
nearest neighbor sites. In later sections we will describe such a model, but we will begin, 
for simplicity, by ignoring the spatial aspects. We do this by imagining that we pluck all 
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the A, B and E from their particular sites and put them into a single large patch, with no 
record of their original spatial locations. Any memory of which individuals were nearest 
neighbors is now lost, and any two individuals picked at random are just as likely to interact 
as any other two similar individuals picked at random. In fact the time evolution of the 
spatial model which would read: pick a site at random, then pick a nearest neighbor of this 
site and implement the interaction rule for the two chosen individuals, now reads: pick two 
individuals from the patch and implement the reaction rule for the two chosen individuals. 

There are other, slightly different, ways of arriving at the above picture. For example, 
instead of dividing up the area of land into sites containing either one individual or no 
individuals, it could be divided up into a number of smaller patches, each of which contains 
several individuals. In this way of thinking, each small patch contains several A, B and 
E types, which interact with each other in the same way as for the large patch described 
above, and which move by exchange interactions with neighboring patches. While this 
defines a slightly different spatial model, the non-spatial version of the model is the same: 
all the individuals from the various patches are collected together into a single large patch. 
Moreover, even the spatial version of the model is in some sense a "coarse-grained" version 
of the original model — several sites in that model when viewed on a coarser scale can be 
reinterpreted as a site in the latter model. These features will be explored in more detail 
when we discuss the spatial aspects of these models in Section |Vl For the remainder of this 
section we will consider only the non-spatial model. 

Suppose to begin with we consider the simpler case of a single species, that is, a patch 
containing only A and E individuals. We shall postulate that the population dynamics of 
the system can be essentially described by three processes: birth, death, and competition. 
The first and third processes will involve two individuals: AE — > AA (birth) and AA AE 
(competition), but the second process involves only one individual: A ^ E (death). These 
seem natural choices since, while death can be modeled as constant, independent of the 
density of individuals, the reduction in the numbers of A due to competition and the growth 
in numbers of A due to births will be density dependent. In other words, there will be a 
tendency for AA to go to AE because of overcrowding, and for AE to go to AA due to the 
presence of resources (space) to sustain a new individual. 

The time evolution of the model can now be described. At each time step we sample the 
patch. On a fraction yU of these occasions we randomly choose 2 individuals and allow them 
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to interact and for (1 — /i) of the draws we choose only one individual randomly. If in the 
former case we draw two E "individuals" or in the latter case one E "individual" , we simply 
put them back into the patch. For all other choices, an interaction may occur leading to the 
replacement of a different set of individuals to those drawn. For each of these processes we 
will introduce rate constants b, c and d as follows: 

AA^AE, AE-^AA, A^E. (1) 

We now only need to know the probabilities of drawing various combinations from the patch. 
Simple combinatorics gives 

n (n — 1) 



Probability of picking AA = /i 
Probability of picking AE = 2/i 



N N-1' 
n [N — n) 
N N-1 ' 



TL 

Probability of picking A = [1 — fi) — , (2) 

where the factor of 2 in the second term comes from the fact that the choices AE and EA are 

identical. These results enable us to write down expressions for the transition probability, 

per unit time step, of the system of individuals going from a state with n A individuals to 

a state with n' A individuals. We denote this quantity by T{n'\n). Since only transitions 

from n to n ± 1 may take place during one time step, the only non-zero T{n'\n) are 

n (n — 1) \ , 

r(n-l|n) = /ic— +{l-fi)d—, 

T{n + l\n) = 2f,b-^-j;^. (3) 

The process defined by (jS)) is a one-step Markov process and so we can immediately write 
down a master equation describing how the probability of having n individuals present in 
the patch, P{n,t), changes with time S^l- The rate of change of this quantity with time 
is simply the sum of transitions from the states with n + 1 and n — 1 A individuals to the 
state with n A individuals, minus the sum of transitions from the state with A individuals 
to the state with n + 1 and n — 1 A individuals: 
dP{n,t) 



dt 



= T{n\n + l)P{n + 1, t) + T{n\n - l)P{n - 1, t) 
- T{n - l|n)P(n, t) - T{n + l\n)P{n, t) . (4) 



This set of coupled equations has to be solved subject to an initial condition, typically 
P{n, 0) = 5„,no! that is, a condition stating that there are known to be Uq individuals in the 



patch at t = 0. Care should also be taken with the boundary values n = and n = N, 
since not all of the transitions are present in these cases. From (jS)) we see that T(— 1|0) and 
T{N + 1|A^) are formally zero. So as long as we define T(0| — 1) = T(A^|A^ + 1) = 0, thus 
we may use the general form (0)) even for n = and n = N . 

The master equation (0)) gives a complete description of the time evolution of the non- 
spatial model. In the next section we will discuss the model predictions in more detail. 
Here we simply wish to make contact with the mean field (i.e. the deterministic) version 
of the model, obtained by taking the N oo limit. This is most easily accomplished by 
multiplying the master equation (j3)) by n and summing over all values of n. By shifting the 
variable in two of the sums on n by +1 and —1, the following rate equation is obtained: 

-^ = Y^T{n + l|n)P(n, t) - ^ T(n - l|n)P(n, t) , (5) 

n=0 n=0 

where angle-brackets signify averages over the possible states of the system. Defining 

AT- 1' N -V N ' ^ ' 

and using Q, ® becomes 

d (n) l/''^/ ~ / n f n l\\ 1 / \ 



So far no approximation has been made in the derivation of ((2)). However we now take 
the limit ^ oo. In addition to eliminating the factor in the second term on the right 
hand-side of ((Tj), it allows us to replace {v?) by {n)"^. This gives 

^ = 2h(t) (1 - 0) - - d(t) , where <P = ^ ■ (8) 

Here (j){t) is the density of individuals in a given area. It is more conventional to write this 
in the form 

^ = NA{r- aNA) , where iV^(t) = {n) = N<f){t) , (9) 



where 



r = 2b-d, a=^. (10) 



Eq. Q is the mean field equation of the model and is the familiar logistic equation, usually 
written down as a phenomenological description of the population growth of a single species 
with intraspecific competition. Here it is derived as the N oo limit of our stochastic model 
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and thus provides a reasonable description of our system when the potential size (number 
of A plus number of E types) of the system is relatively large. Of course, this limit is 
purely formal. In practice what we mean that if is of the order 10^, for instance, then this 
approximation is good if we are only interested in accuracies of up to 0.01% (if the next order 
corrections are of order 1/A^) or 1% (if the next order corrections are of order 1/^/N). This 
approximation obviously cannot describe chance extinctions, which occur when n is small, 
nor does it predict a mean time to extinction for the A population. In regimes where these 
effects are important, it provides a poor description of the system, but this will inevitably 
be true of any purely deterministic description. 

The necessity of introducing the "empty site" individuals E should be clear from the 
above derivation. In order to be able to define a population density which changes with 
time we need a null population which can be displaced if the A species is successful and 
increase if the A population falters. It is also very natural to have a ceiling on the growth 
of A individuals {Na < A^) representing a limit on the available resources. If no E's were 
introduced in the two species case, the two population sizes would not be independent and 
would simply add up to A^. The A population would obey (jH)) and Nb = N — Na- This 
is clear if we simply imagine repeating the above derivation for a two species system by 
replacing E by B. Although the interaction rules would be altered by this replacement, 
it would still be the case that A^^ + A'^ = A^. So it is vital to have "empty space" for 
individuals to exploit in order to obtain a realistic population dynamics. 

We have discussed the single species case in some detail since the construction of the 
master equation in the two species case (and, in fact, in the S-species case for arbitrary S) 
follows similar lines. We still draw 2 individuals fi of the time and 1 individual (1 — /i) of 
the time, but the processes and their rate constants are now: 

AA^AE, AB^AE, BA^BE, A^E 

BB^BE, AE ^ AA, BE ^ BB, B^E. (11) 

The rate constants Caa, a = 1, 2, represent intra-specific competition and Cajj a f3, inter- 
specific competition. Analogous probabilities to of choosing particular combinations of 
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A, B and E are 

n in — 1) , n (N — n — m) , , n 

m (m — 1) m (A^ — n — m) , , m 

Transition probabilities now have initial and final states specified by two integers. The 
transition probability per unit time from the state {n, m) to the state (ra', m') will be denoted 
by T{n' ,m'\n,m). The non-zero transition probabilities are 

Ti(tI — Tl TTITI 

T{n - 1, m|n, m) = fj, Cn — — ■ + (1 - /x) rfi — + 2/^ cu 



N{N-1) ^ ^iv ' ^ '^A^(iV-l)' 

I m(m — 1) \ , ^ ''^'^ 

r(n, m - l|n, m) = /i C22 ^77777 -r + (1 - /i) 02 — + 2/i C21 



+ 1, m|n, m) = 2/i6i 



r2(A^ — n — m) 



N{N-1) ' 

m(A^ — n — m) . . 

T(n,m+l|n,m) = 2/162 _ — • (1^) 

The master equation is an obvious generalization of (jH): 

— ^ — = T(n, m|n + 1, m)P(r2 + 1, m, t) + T(n, m|n — 1, m)P(n — 1, m, t) 

+ T{n, m\n, m + l)P(n, m + 1, t) + T(?7., m|?7,, m — l)P(n, m — 1, t) 

— {T(n — 1, m|n, m) + T(n + 1, m|n, m) + T(n, m — l|n, m) 

+ T(n, m + l|n, m)} P(n, m, t) . (14) 

This equation simply expresses the increase in P(n, m.t) due to the four possible transitions 
into the state (n, m) described by T(n, m|r2 ± 1, m ± 1), and the decrease due to transitions 
out of this state described by T(n ± 1, m ± l|n, m). The boundary and initial conditions are 
obvious analogs of those in the one-species case. The generalizations of (0) are quite simple, 
since none of the transition probabilities involving changes only in m enter into the equation 
for d{n)/dt, and none of the transition probabilities involving changes only in n enter into 
the equation for d{m)/dt: 

din) ^ ^ 

= T{n + l,m\n,m)P{n,m,t) — T{n — l,m\n,m)P{n,m,t), 

n,m,=0 n,m=0 



d{m) 



N N 



T{n,m + l\n,m)P{n,m,t) — T{n,m — l\n,m)P{n,m,t) . (15) 



dt 

n,m=0 n,m=0 
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We can now substitute the forms for p3p into this equation and take the mean field hmit 
oo. This allows us to factor {mn) into {m){n), as well as {n"^) into (n)^ and (m^) 
into (m)^ as before. The final result for the rate of change of population densities are the 
competition equations 

dNA 



dt 
dNB 



Na (ri - aiiNA - auNs) where Na = {n), 
Nb {^2 - a2iNA - a22NB) and where Nb = {m) , (16) 



dt 

familiar from population biology textbooks. The parameters in ()l(j|) are related to the 
parameters of the stochastic model by 

2flba {1 - fl)da fl{2ba + Caa) 2/x(6a + Ca/3) 



N-1 ^ . N{N-1) ' iV(iV-l) 

As we will see in section essentially the same kind of reasoning as that given above 
can be applied in the spatial version of the model. Before discussing this, however, we will 
investigate the large N limit of (jH) and (fT^ a little more carefully, obtaining corrections to 
mean field theory and comparing these to simulations. 

III. BEYOND MEAN FIELD THEORY FOR THE NON-SPATIAL MODEL 

In the last section we gave arguments to show that the mean field versions of the stochastic 
models we had introduced were indeed the deterministic models conventionally used to 
describe these systems. In this section we will apply an elegant method due to van Kampen 
37| which not only allows us to obtain these results in a more systematic way, but also gives 
a method of finding stochastic corrections to this deterministic result for large N. We will 
only work to next-to-leading order in this paper. This will give a Gaussian broadening to 
P[n,t), or P{n,m,t), with the peak of the distribution moving according to the relevant 
deterministic equation. We will then compare these results with numerical simulations of 
the full stochastic process. The large expansion is very clearly explained by van Kampen 
in his book [s^j, so we will content ourselves with giving a brief outline of the method as 
applied to the one species case. The two species calculation follows very similar lines. 

We saw in the last section that, in the limit N ^ oo, the variable n became deterministic 
and equal to N(j){t). In this limit the function P{n,t) will be a delta-function. For large, 
but finite A^, we would expect P{n,t) to have a finite width of order N.N~^/'^ = N^/'^. Now 
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n is once again a stochastic variable, and it is natural to bring out the large structure of 
the theory by transforming to a new stochastic variable ^ by writing 

n = N(j){t) + N^/'^^. (18) 

We will not need to assume that the function (j){t) satisfies any particular differential equa- 
tion; if we simply choose it to follow the peak of the distribution as it evolves in time, then 
the equation it satisfies will emerge. A new probability distribution function 11 is defined 
by P{n,t) = Il{^,t), which implies that 

dt dt di ^ ' 

When using this formalism it is useful to rewrite the master equation (j3]) using step operators 
which act on an arbitrary function of n according to £f{n) = f{n + 1) and S~^f{n) = 
f{n— 1). This gives 

dP{n,t) 



dt 



{S - 1) [T{n - 1 \n)P{n, t)] + {S~^ - l) [T{n + 1 \n)P{n, t)] . (20) 



This form of the master equation is useful because the step operators have a simple expansion 
involving powers of the operator N~^^'^d/d^, which simplifies the identification of differing 
orders in N'^^"^ . We shall assume the initial condition on the equation to be P(n, 0) = 5n,no- 
Applying the method, and identifying powers of N^/^, yields the macroscopic equation 

^ = «i,o(0) (21) 
to leading order and a Fokker-Planck equation 



^ = -<o(^)^CT + -«,o(^)^ 



"Io(</')t^ + -«2,o(</')i^ (22) 



describing a linear stochastic process to next order. Here the functions ai^o(0) and 02,0 (0) 
(we have used van Kampen's notation) are given by 

ai,o(0) = 260(1 -0) -0(^+50) 

o;2,o(0) = 260(1 - 0) + 0(J + c0) . (23) 



Since the Fokker-Planck equation (j22|) describes a linear process, its solution is a Gaussian. 
This means that the probability distribution n(^,t) is completely specified by the first two 
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moments (^)t and Multiplying ()22|) by ^ and and integrating over all ^ one finds 

dtiOt = «'i,o(0)(Ot 

dt{e)t = 2<o(0)(ai + «2,o(0). (24) 

The procedure is to solve (PT|) and obtain as a function of t. This function is then 

substituted into and these equations solved for (^)t and {^^)t- In the case of the first 
moment this may be performed quite generally to give 

(Ot = (Oo exp |^*rf™l,o(0W)} (25) 



Choosing our initial condition to be 



m 



no 

N 



(26) 



the initial fluctuations vanish, and (.^)o = 0. Thus {C,)t = for all t. 

This summarizes the method. We therefore start by solving ()21|) . subject to (f^Bj) . Defining 



p = 2b — d ; a = 2b + c, 
for convenience, the general solution for p 7^ is 



0(t) 



P 



where the constant A is determined by the initial condition 

P 



A = a- 



m' 



(0(0) ^ 0) 



(27) 



(28) 



(29) 



If 0(0) = then 0(t) = for all t. If p = a degenerate form of the solution exists and is 
given by 



0(0) 



(p = 0) . 



(30) 



1 + (70(O)t' 

This solution can now be substituted into the equation for (^^)t given in ()24j] . An inte- 
grating factor for this equation is e^''Y0^(t), which yields upon use of the initial condition 
(e')o = when p ^ 0, 

1 



\a - Ae-P*] 



- a A 



4b^ + 106(c + d) + cd e-P' [I - e-P'] 
2A^p 4b^ + 4b{c + d) + cd te-^P' - A^{2b + d)e-'^P' [l - e"^*] } (31) 
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Of course, care has to be taken when applying this approximation. If the distribution has 
significant interaction with the boundaries ai n = or n = N, the Gaussian approximation 
will break down. So, for example, if p > 0, the peak of the distribution will move from uq 
eventually coming to rest at Np/a. While this is happening the probability distribution 
broadens, eventually reaching its stationary value at 

' (26 + £)2 ^ ' 

On the other hand, if p < 0, the peak of the distribution will eventually tend to zero, and 
so the Gaussian approximation will break down at some finite time. 

The case of two species follows in an exactly analogous manner: one writes n = N(j){t) + 
N^^"^^ and m = Nip(t) + N^^'^rj and defines a new probability distribution 11 via P(n, m, t) = 
r], t). The macroscopic equations obeyed by (f){t) and ipit) are again found to be identical 
to those given in section |nj and 11 is again found to satisfy a Fokker-Planck equation 
describing a linear process, and is therefore a multivariate Gaussian distribution. Details 
are given in Appendix A, where one sees that analytic expressions for the analogs of 
cannot be obtained, in part at least, because the macroscopic equations cannot be solved in 
closed form. However, we have solved them numerically, and we now go on to compare the 
large results in both the one species and two species cases with simulations of the original 
stochastic process. 

IV. SIMULATIONS OF NON-SPATIAL MODELS 

In order to better understand the range of validity of mean-field theory and its Gaussian 
corrections, we have performed numerical simulations of the non-spatial model described 
above. The simulation of stochastic models of population dynamics has progressed in line 
with the availability of high-speed computers over the last two decades. During the beginning 
of this period, books on the stochastic dynamics of fluctuations in biological systems only 

Ssj l- about a decade ago they had assumed a more central 



mentioned simulations fieetingly 



role [39|, and are now regarded as essential in the understanding of these systems |40|, 141 1. 
We refer the reader to these last two references for more details on how these simulations 
are carried out in practice. 

In our numerical algorithm an ensemble of patches is iterated forward in time. To achieve 
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reasonable statistics we generally take the size of the ensemble to be several thousand real- 
izations. In each small time step a single individual or a pair are selected from each patch 
in turn and with the appropriate probabilities transitions are made and the new individuals 
replaced. 

We periodically average over the ensemble to measure both the mean densities of in- 
dividuals and the variance in the densities. Concurrent with the stochastic simulation we 
also integrate forward the mean field equations and the equations for the Gaussian vari- 
ances (Eqs. (j23) and ()A5|) - ()A9|) ) to allow a direct comparison. Forward integration of the 
differential equations is performed using a second order Runge Kutta scheme. 

We will present some examples of our numerical work which illustrate the main effects 
that we have found. We have not performed an exhaustive numerical analysis due to the 
large parameter space of the competition models. The examples we give here are fairly 
typical of a wide range of parameter space, and are chosen to illustrate a variety of effects. 
Our results may be summarized by the statement that so long as the size of the patch is large 
and the system is not close to extinction (such that discrete effects play a strong role) then 
the mean field equations, and the large correction yield a remarkably accurate description 
of the system dynamics. However, there is a fairly sharp transition signaling the failure of 
mean field theory as the patch size is reduced below a critical value. For smaller patches, 
the quantitative precision of the mean field theory fails badly. This inaccuracy gives way 
to qualitative errors if one runs the system close to extinction. In this case the probability 
distribution of the populations is poorly approximated by a Gaussian function and one is 
compelled to abandon mean field theory and its Gaussian corrections. It is significant that 
the critical patch size depends sensitively on whether the patch contains one or two species, 
and whether there is interspecific competition. To clarify these statements we now lead the 
reader through some illustrative examples. An exhaustive list of parameter values is given 
in Tables 1 and 2. 

In Figure 1 we give an example of a moderately large patch (A^ = 100) containing a 
single species. The mean density soon settles down to a (quasi) steady-state value as does 
the variance in the population density. Mean field theory and the large A^ corrections give 
very good agreement. In Figure 2 we show results for an identical situation but with the 
patch size reduced from 100 to 10. In this case the quasi-steady-state is meaningless since 
extinction events are frequent and the mean density (measured over the ensemble of patches) 
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steadily decays to zero j38|, |42]. Note, that the variance also decays to zero, since as time 
proceeds more and more realizations go extinct and the probability distribution of population 
densities is dominated by a delta function peak at zero. This illustrates the effect of discrete 
individuals for small systems. 

We now consider competition between two species in a single patch. We only con- 
sider large patches, where naively, discrete individual effects may be neglected. In Fig- 
ure 3 we study the simplest case in which the A and B individuals have identical birth, 
death, and intraspecific competition rates. There is no explicit interspecific competition 
(i.e. Ci2 = C21 = 0), although the finiteness of the patch leads to indirect competition be- 
tween all individuals (see Eqs. and (|17p). The patch is taken to be very large with a 
capacity of 400 individuals. We find satisfactory agreement between simulations and the 
mean field theory and its large-A^ corrections. It is somewhat surprising that on reducing 
the patch size from 400 to 200 (Figure 4) the variance in the system increases steadily, 
far exceeding the large corrections. Returning to the large 400 capacity patch and now 
introducing a small amount of interspecific competition (Figure 5) we again find that the 
large corrections fail to capture the growing fluctuations in the system. We conclude from 
this, and other simulations, that mean field theory and its Gaussian corrections can work 
well, but only for patches above a critical size. This critical size is itself strongly dependent 
on the number of species, various growth parameters, and the presence or not of explicit 
interspecific competition. 



V. SPATIAL MODELS 

We have already discussed the spatial versions of the model in Section IHI In one version 
of the model, the area under consideration is divided into a large number of patches, each 
containing a small number of individuals, which are then identified with the sites of a 
regular two-dimensional lattice (usually a square lattice). Competition takes place between 
individuals of a particular patch, and the birth rate is similarly only dependent on the 
population density of the parental patch, but individuals are allowed to migrate to nearest 
neighbor patches, if space is available (that is, if an empty space, E, exists at the neighboring 
site). In terms of this picture, the lattice consists of an array of patches, which interact 
through migration of individuals from one patch to a nearest neighbor patch. In the other 
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version of the model introduced in section m each patch contains only one individual, thus 
the sites of the lattice represent individuals rather than patches. In this case competition 
and birth processes, as well as migration, depend on the occupancy of nearest-neighbor sites. 
In both versions of the model, the death rate is constant. 

It is clear that many other variants are possible. In general, the first model will be 
more applicable to situations where individuals move on length scales which are much larger 
than the communities they live in, and the second model more applicable when all these 
processes occur on scales which are of the same order. However, we will see that in the 
exploration of the nature of the mean-field limit and the importance of stochastic effects, 
which is what interests us in this paper, these differences play a secondary role. In the next 
two sections we will follow the same program as was carried through for the non-spatial case: 
describing the stochastic process which defines the model, obtaining the mean- field limit, and 
finally comparing the mean field equations with simulations of the full model. Much of the 
mathematical detail will be relegated to Appendix B; unfortunately while many of the ideas 
are simple generalizations of those introduced earlier, the mathematical notation becomes 
(of necessity) rather complex and detracts from the points which we wish to emphasize. 

The simplest spatial process to describe is that of a single species in the first version of 
the model. The processes in this case may be broken down into 3 classes: 

(i) For a fraction qi of the events we randomly pick a site i and then randomly draw two 
individuals from within the patch at that site. If two E individuals are drawn, they 
are simply replaced, otherwise the following interactions may occur (c.f. Eqn. (^) 

AAi AE,, AE, AAi. (33) 

(ii) For a fraction q2 of the events we randomly pick a site i and then randomly pick 
another site j which is a nearest neighbor of i. One individual is drawn from the patch 
at i and another from the patch at j. If these two individuals are of the same type 
(two A's or two -E"s) no action is taken, otherwise migration with a rate constant m 
may occur: 

AEj ^ EiAji EiA, ^ AE,. (34) 

(iii) For a fraction 1 — gi — g2 of the events we randomly pick a site i and then randomly 
draw a individual from within the patch at that site. If an E individual is drawn no 
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action is taken, otherwise death may occur at a constant rate d: 



(35) 



The probabihties of choosing these various processes are in the case of (i) and (iii) simply 
modifications of Q. The modifications that are required are that all n should be written 
as rii to denote the number of A's in the patch at site i, /i should be replaced by qi, (1 — /i) 
by (1 — gi — ^2) and all terms multiplied by fi^^, where Q is the number of sites in the 
lattice. We assume that the number of individuals in each patch is the same for all sites 
and is denoted by A^. If the sites i and j have already been chosen, the probabilities for the 
processes (ii) are: 

Ui {N - Uj) 



Probability of picking AiEj = q2 



N N 



Probability of picking EiAj = — — • (36) 

N N 

The transition probabilities, master equation, and the derivation of the population-level 
equation corresponding to this model, are discussed in Appendix B. The lattice version of 
this latter equation is given by ()B8|) . On taking the continuum limit, defined by ()B9|) . it 
becomes 

^ = mV^^ + 2h(P (1 - 0) - c(p^ - d(P , (37) 

OT 

where (/)(x, r) is the continuum version of {ni{T))/N in the limit N 00 and where r is 
a rescaled time. This equation is exactly the result (jHI) obtained in the non-spatial case, 
but with the addition of a drift term. We can write (|37p in a more standard form by 
defining a diffusion constant Da = rh and making the identification ()1U|). This gives 

= DaV^Na + iVA (r - uNa) , (38) 

where, as before, Na = N(j). 

The discussion of the second version of the spatial model follows similar lines. Now there 
is only one individual per site, and therefore rij can only take on only two values: and 1. 
In addition, birth and competition processes, as well as migration, depend on the occupancy 
of nearest neighbor sites. Therefore there are only 2 classes of processes: 

(i) For a fraction fi of the events we randomly pick a site i and then randomly pick another 
site j which is a nearest neighbor of i. If these sites are both E^s no action is taken, 
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otherwise migration with a rate constant m may occur according to ()34p. or birth and 
competition may take place with rate constants b and c/2 respectively: 

AE, AAj, EiA.j AiAj, AiA, ^ EiA.j, A^A, ^ A^E,. (39) 

The factor of 1/2 has been introduced into the rate constant for competition in order 
be consistent with the non-spatial case and the first version of the spatial case. 

(ii) For a fraction 1 — /i of the events we randomly pick a site i. If the site contains an E 
individual, no action is taken. Otherwise death may occur at a constant rate d given 

by dsni). 

Since there is only one individual per each site, the probabilities of picking AiEj or EiA^ 
are simply pUj) with = 1 and q2 replaced by /i. Similarly, the probability of picking Ai 
is as in the first version, but again with = 1 and with (1 — (?! — q'2) replaced by (1 — /i). 
The only new feature is: 

Probability of picking AiA^ = finiUj . (40) 

Just as before, we have assumed that the sites i and j were already chosen, so that the above 
probabilities only represent the choices of types of individuals at these chosen sites/patches 
and also the choice of the number of individuals in an event (one or two). In the first model, 
we denoted the number of sites in the lattice by Q. This was independent of A^, the number 
of individuals in a patch. It was this latter quantity that we allowed to become infinitely 
large, in order to deduce the population-level description. In this second version, there is 
only one individual per site, and so it is the number of lattice sites (now denoted by A^) 
which we take to be infinitely large. More details of this approach are given in Appendix B, 
where it is shown that, in the large A^ limit, the population-level description is again given 
by ()37p — albeit with slightly different definitions of the parameters. This is not a surprise; 
we would expect there to be a large number of IBMs which differ in detail, but which have 
the correct qualitative features, and give the same population-level description. 

The partial differential equation (jHH|l is simply the ordinary differential equation for the 
non-spatial case (jH)), but with a term V^A^i added. So the corresponding spatial description 
is indeed obtained by using the simplest prescription. However, this will not turn out to 
be the case when more than one species are present. It is this scenario which is of most 
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interest to us in this paper; we have described the one species case in some detail, largely 
because it is technically simpler, and therefore the crucial steps in the argument clearer. The 
many species case may differ at the population-level, but the the setting up of the IBMs 
and the derivation of the population-level equations are a straightforward generalization of 
the one-species case. 

Let us once again begin with the first version of the model where birth, competition, 
and death processes are purely local (they take place in a single patch at a specific site on 
the lattice) and only the process of migration involves nearest-neighbor patches. All of the 
transitions are variants of those in models previously considered in this paper. Specifically 
the 3 classes of processes are: 

(i) For a fraction qi of the events we randomly pick a site i and then randomly draw 
two individuals from within the patch at that site. If two E individuals are drawn, 
they are simply replaced, otherwise the interactions are given by the 6 two individual 
interactions in with a site index i added on to the A, B and E individuals. 

(ii) For a fraction q2 of the events we randomly pick a site i and then randomly pick 
another site j which is a nearest neighbor of i. One individual is drawn from the 
patch at i and another from the patch at j. If neither of these two individuals are -E"s 
(no space) or both are £"s (no migration possible), then no action is taken, otherwise 
migration with rate constants mi or m2 may occur: 

AiEj EiAj, EiAj AiEj, BiEj EiBj, EiBj BiEj. (41) 

(iii) For a fraction 1 — gi — g2 of the events we randomly pick a site i and then randomly 
draw a individual from within the patch at that site. If an E individual is drawn no 
action is taken, otherwise death may occur at constant rates di or d2'- 

A, ^ E,, B, ^ E,. (42) 

The probabilities of choosing these various processes are in the case of (i) and (iii) simply 
modifications of ()12j) . The modifications that are required are exactly those we described in 
the similar version of the one species case: the n's and m's should be written as rii and 
respectively, /i should be replaced by qi, (1 — /i) by (1 — gi — ^2) and all terms multiplied 
by Q^^. The migration of A's and 5's are independent of each other, and so are described 
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in exactly the same way as for single species. Details are given in Appendix B where it is 
shown that, after the continuum limit has been taken, the equations for 

0(x,r) ^ hm and ^(x,r) ^ lim Mill (43) 

are 

+ 26i0 {l-(P-ij)- gn(/.2 + 2£i2(/'V^ - di^ , (44) 

and 

+ 262^ (1 - - ^A) - £22^' + 2c2i^/'(/' - • (45) 

The second version of the two-species model has rij = 0, 1 and rrii = 0, 1, with birth, 
competition, and migration depending on nearest-neighbor occupancies. The 2 classes of 
processes are 

(i) For a fraction fi of the events we randomly pick a site i and then randomly pick another 
site j which is a nearest neighbor of i. If these sites are both £"s no action is taken, 
otherwise migration may occur according to ()41|). birth according to 

AiEj AiAj, EiAj AiAj, BiEj — ^ BiBj, EiBj — ^ BiBj, (46) 

and competition according to 

AiAj > AiEj, BiAj >■ BiEj, AiBj — > AiEj, BiBj — > BiEj, 



AiAj EiAj, AiBj EiBj, BiAj EiAj, BiBj EiBj. (47) 

(ii) For a fraction 1 — /i of the events we randomly pick a site i. If the site contains an E 
individual, no action is taken. Otherwise death may occur according to 

The probabilities of picking two individuals, one of which is an E, are the same as in the 
first version of the model, but with = 1 and q2 replaced by fi. The probabilities of picking 
a single individual are similarly related to those found in the first version. The probabilities 
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associated with picking AiAj is given by ()4U|). BiBj by fimimj and AiBj by finirrij. In 
Appendix B we describe how, in the hmit where the number of lattice sites N becomes 
infinitely large, the continuum versions of (rii) and {nii) again satisfy equations and 
(gSl). 

So, in summary, both versions of the IBMs we have discussed in this paper give rise to 
the same population-level equations. This is true whether there is only a single species in 
the system, or whether two species are present. In the one species case this equation is given 
in standard form by ()38p. To write the two species equations ()44j) and in standard form 
we make the identification (|T7|l and introduce diffusion constants 

D^ = mi, DB = m2, ^1 = ^, ^2 = ^- (48) 

This gives 

DaV^Na + Di {NaV^Nb - NbV^Na) + Na {n - GuNa - a^^NB) , (49) 



dNA 

dr 



o TV r 

= DbV^Nb + D2 {NbV^Na - NaV^Nb) + Nb (ra - a2iNA - asaiVs) , (50) 

where Na = N(j) and Nb = Nifj. Unlike ()38|1 . these are not the standard equations found in 
population biology textbooks. 

The additional terms which appear in pUj) and (jM^j) . but not in the standard equations, 
are antisymmetric in A^^ and Nb and involve derivatives, and so do not appear in non-spatial 
models or spatial models with only one species. Their structure is dictated by the way that 
migration is modeled at the individual level. Since their occurrence is generic, they will also 
appear in spatial models derived from IBMs having three or more species. Although these 
terms have not to our knowledge been discussed in the context of ecological models, they are 



well-known in the context of interspecies diffusion 
in quantum field theory 



4J| in physics, and they also appear 



VI. SIMULATIONS OF SPATIAL MODELS 



The simulations of spatial competition was performed in a directly analogous fashion to 
the non-spatial model described in Section HVl We confined ourselves to one spatial dimen- 
sion for simplicity. In addition, we only simulated the version of the model where a patch of 
size A^ was placed at each site of the one-dimensional lattice of size L. Within each patch 
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the usual dynamics of competition is played out. Furthermore, in each small time step there 
is a small probability of dispersal of individuals from a given patch to the two neighboring 
patches, just as encoded in the master equation. We generally started with initial conditions 
in which species A and B were spatially separated and then proceed to intermix and com- 
pete as individuals diffuse from patch to patch. The mean field equations (|49p and (j5U|) are 
again integrated forwards in time using second-order Runge-Kutta methods. An exhaustive 
list of parameter values is given in Table 3. 

In a spatial system such as this, extinction is much less of a problem since should a patch 
become empty it will soon be restocked from neighboring patches. Despite the weakened 
effect of discreteness in small patches, we still find that the behavior of the spatial systems 
differs significantly from mean field theory when patch sizes are below a critical value (of 
approximate value 50 for the results presented here). As before, we have chosen to present 
two typical scenarios. 

In Figure 6 we show the early and late time behavior for a system in which initially the A 
individuals occupy the left half of the system, and the B individuals occupy the right half. 
In the ensuing dynamics, the A and B individuals have identical mobilities, growth rates, 
and intraspecific competition parameters. However, the A individuals are disadvantaged 
by having a slightly higher death rate than the -B's. This is balanced by giving the A's an 
interspecific competitive advantage over the 5's. On varying the strengths of these balancing 
forces it is possible to obtain invasion of A's from left to right or invasion of -B's from right 
to left. We have chosen an example of the latter. It is seen that mean field theory does an 
excellent job in predicting the long time dynamics of the system. In this figure the patch 
size is rather large with a capacity of 100. 

In Figure 7 we repeat the exact simulation as before but simply reduce the patch size 
from 100 to 10. In this case the A individuals are severely affected by discrete extinction 
events and their density is in poor agreement with mean field theory. Interestingly, the 
denser B individuals are fairly well described by mean field theory throughout the range. 

We also studied an alternative balance of effects as follows. In Figure 8 we show a 
situation in which the death rates for the two species are the same, but now we reduce 
the amount of intraspecific competition among the B individuals. Again, the A's have an 
interspecific competitive advantage over the -B's. In this case the A's invade the -B's. Here 
the patch size has the intermediate value of 50 individuals. We see that mean field theory 
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performs relatively well. 

On reducing the patch size for this particular scenario, from 50 to 10 (Figure 9) we see 
the failure of mean field theory (which predicts invasion from left to right). The enhanced 
fluctuations in the smaller patches lead to a quasi-dynamical balance in the interfacial region 
between A's and -B's. In this region the A's are beset by fluctuation- induced extinction events 
and this makes them too weak to invade the i?'s in the usual manner of a Fisher wave. 
Instead, over longer scales than shown in the figure, the density of A's slowly permeates the 
5-rich region in a "creeping" motion. 

VII. CONCLUSIONS 

There are many ways to formulate population dynamics. Popular descriptions tend to 
be either deterministic (mean field) equations, or individual based algorithms designed for 
implementation on a computer. The extreme difference in these two approaches has led to 
difficulties in directly comparing results. Disparities may be due to fundamental deficiencies 
in one or both of the techniques, or else be attributable to "renormalization" of various 
parameters. In this paper we have attempted to bridge the gap between mean field models 
and individual-based models. We have described a very general framework with which to 
formulate population dynamics using the language of "patches" to create a concrete picture 
of the stochastic process. The size of the patch is the central parameter. Mean field theory 
is recovered on taking the patch size to infinity, while discrete stochastic effects become 
prominent for small patches containing a few individuals. Again, we emphasize that in our 
usage "mean field" refers to the approximation in which cross-correlations between stochastic 
variables is neglected, but still allows for an explicitly spatial description. 

From a biological perspective, a patch can be thought of as a (small) spatial region within 
which interactions between individuals occur. It is assumed that movement within this scale 
is not biologically significant. In our spatial patch model, movement of an individual between 
patches is biologically significant since that individual will now have interactions with a new 
set of individuals in a neighboring region. For systems in which interactions (not involving 
movement per se) occur over larger scales, it will be necessary to include additional inter- 
patch processes. 

We have studied both non-spatial and spatial models. The non-spatial case corresponds 



25 



to a single patch containing a number of individuals of both species. We have derived the 
corresponding mean field theory and its first order corrections (i.e. Gaussian fluctuations 
about the deterministic predictions). Generally, so long as the patch size is above a critical 
value (which tends to be of the order of 100 in the examples shown here), and the birth and 
death rates are such that a sizable quasi-steady-state population is possible, then the mean- 
field theory and its corrections give a satisfactory description of the system. For smaller 
patches, or for situations in which there is a non-negligible probability of extinction, it is 
crucial to account for the discrete nature of the individuals. The population dynamics is 
inherently stochastic and one must dispense with a deterministic description. By tuning 
the patch size we have seen that the transition from a mean-field like to a stochastic regime 
is rather sharp and dependent on the existence of inter-specific interactions (in this case, 
competition). 

The same general picture holds for the spatially explicit models. We have discussed two 
types of spatial patch models. In one, at each spatial site there is a "micro-patch" which 
may hold at most one individual. Movement and competition occurs between patches. In 
the other, each lattice site is a patch of tunable size and competition occurs inside the 
patch. Movement, of course, is still between patches. A careful formulation shows that 
each model has the same spatial mean- field limit. Of particular interest is the emergence 
of novel non-linear diffusion terms, which are only present when two or more species are 
competing for space. These terms are not written down in the standard "intuitively derived" 
continuum equations of spatial competition models. They are especially important in spatial 
regions in which the density of one species is high, while the density of the other is strongly 
spatially varying. This would occur, for instance, in a region of space containing a population 
boundary for one species (due to some environmental barrier) but not for the other. We 
intend to investigate such effects in more detail in a follow-up paper. 

In our investigation of spatial mean field models, we have found them to be more robust 
than in the non-spatial case. This is primarily due to the weakening of local extinction via 
continual rescue effects from neighboring patches. It is still the case however that as the 
patch size is decreased, the quantitative precision of mean-field models suffers, and with 
smaller patches still (we have in mind patches of size 10 or less) new stochastically driven 
qualitative features emerge. An example of this was given (Figure 9) in which an invasion 
process (in mean field theory) was halted due to stochastic weakening of the leading edge of 
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the invading population density. 

In conclusion, we have presented a simple framework with which to discuss fluctuation 
effects in population biology. This framework is based on the use of patch models as concrete 
realizations of stochastic processes. The transition from mean field behavior to fluctuation 
dominated stochastic dynamics is effected by reducing the size of the patch. The critical 
patch size separating these two regimes depends sensitively on the biological interactions 
present. This has been an intensively theoretical work. In future work we intend to apply 
the patch model to a variety of multi-species population dynamics to address the importance 
of fluctuations and validity of mean field theories in a quantitative and controlled manner. 
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APPENDIX A: LARGE ANALYSIS 

This appendix contains the details of the large analysis for the two species case which 
was described in section UTTl in the one species case, 

It is once again useful to write the master equation ()14|) in the form 

1) [T{n — l,m\n,m)P{n,m,t)] 

- l) [T{n + l,m\n,m)P{n,m,t)] 
1) [T{n,m — l\n,m)P{n,m,t)] 

- l) [T{n,m + l\n,m)P{n,m,t)] (Al) 

where the step operators S are defined by their actions on functions of n and m by 
S^^f{n, m, t) = f{n ± 1, m, t) and S^^f{n, m, t) = f{n, m ±l,t). 

Writing n = N(j){t) + N^/'^^ and m = Nijj(t) + N^^'^rj^ van Kampen's method yields the 
macroscopic equations 

^ = «i,o(0,^) ^=A,o(0,^) (A2) 



dP{n,m,t) _ ,r. _ 
dt ~ 

+ (C^ 

+ {£y - 
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to leading order and the linear Fokker-Planck equation, 
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1 -^^n 1 d^u 

2«2,o ^ + 2/52,0 ^ 



(A3) 



to next order. This is a multivariate Fokker-Planck equation, but it is again linear, and so 
its solution is a (multivariate) Gaussian. The a and /3 functions are given by 



ai,o(0, i^) = 26i0(l - - ^) - <^ C1102 + Ji0 + 25 



A,o(0, V^) = 262^(1 - - ^) - C22^2 ^ ^ 26 



-21 1 



"2,0(0, V') = 26i0(l - - V') + 1 cii0^ + M + 25 



-12V 



/32,o(0, ^) = 2fe2V'(l - - ^) + 622^' + d2ip + 25 



-21^ 



(A4) 



Since the solution to the Fokker-Planck equation is a Gaussian, we need once again only 
find the first two moments. They satisfy 
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(A5) 
(A6) 
(A7) 
(A8) 
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We set the initial conditions on the macroscopic equations by asking that 



(A9) 



(AlO) 



This implies .^(0) = and //(O) = 0, and by successive differentiation of the macroscopic 
equations, that all derivatives of (^)t and {r])t at t = are also zero. We therefore take 

iOt = ivh = (All) 

for all t. Since the macroscopic equations with initial conditions ()A10|) cannot be solved in 
closed form neither can the equations for {C,^)t, {v'^)t or {C,ri)t- 
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APPENDIX B: SPATIAL MODELS 



In this appendix we give details of the transition probabihties and the master equations 
for the spatial models discussed in section |3 of the main text. The results are frequently 
fairly straightforward generalizations of those found for the non-spatial model, however there 
are some surprises in store: for example, the non-trivial spatial terms found in the mean 
field theory of the two species model is only found by a careful step-by-step derivation of 
the equation satisfied by d{ni)/dt. 

We begin with the first version of the one-species model. The transition probabilities for 
the processes defined by and are, by analogy with (jH)), 

rri, II N 9ic rii {ui - 1) (1 - gi - q2)d Ui 
J . . . n,- — 1 . . . . . . n,- . . . = , 

^ " ^ ' ^ n N N-1 n N' 

rri ,11 ^ rii {N -m) 

J ... n,- + 1 ...... n,- .. . = . Bl 

^ ' ' ' ^ n N N-i ^ ' 

The only change is the addition of the factor Vt^^, where Vt is the number of sites in the 
lattice, which represents the arbitrary choice of the lattice site. In the case where the 
transition probabilities involve two neighboring patches, which is the process described by 
the corresponding quantities are 

_ q2'm Ui {N - rij) 



T(. . . rii — 1, rij + 1 . . . \ . . . Ui, rij . . 



zn N N 



T{...n^ + l,nj-l...\...n„nj...) = j;^ , (B2) 

where z is the coordination number of the lattice (the number of nearest neighbors of any 
given site), and represents the choice of the nearest neighbor j, once i has been chosen. 
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The master equation for this process therefore reads 

i T{. . .ni,nj . . .\. . .ni + l,nj -1 . . .)P{. . .rii + l,nj - 1 . . . ,t) 

i jei 

+ r(. . .riijUj . . . \ . . .rii — l,nj + 1 . . ■)P{. . . rij — 1, + 1 . . . , t) } 



+ J2 { T{...ni...\...ni + l...)Pi...n, + l...,t) 

i 

+ T{...n,...\...ni-l...)P{...ni-l...,t) } 

- ^ ^ { T{. . .Hi - l,nj + 1 . . . \ . . .ni,nj . . .)P{. ..ni,nj..., t) 

+ T{. . .rii + l,nj - 1 . . . \ . . .ni,nj . . .)P{. . .ni,nj . . . ,t) } 

-J2 { T{...n,-l...\...n,...)P{...ni...,t) 

i 

+ r(...n, + l...|...ni...)P(...n,...,t) }. (B3) 

Although this looks rather complicated, it is a straightforward generalization of (0)). In 
an effort to keep it as simple as possible, only the number of individuals at sites where 
changes occur [i or j) have been explicitly shown on the right-hand side of the equation. 
The notation j G z denotes a sum over all sites j which are nearest-neighbors of i. On the 
left-hand side of the equation, n denotes the number of individuals in the set of all patches: 
n = (rii, ...,ni,...,nj,.. .)• 

To obtain the rate equation, we substitute ()B3|) into 



_J2n/^f'\ (B4) 



dt ^ dt 

{n} 



Defining new quantities 



QiC I qib 7 {l-qi-q2)d . 

d = 77 , m = — — , (B5) 



N-V N-V N ' N 

as in (jUj), and introducing a rescaled time r = t/Q, the following equation is found: 
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where we have used the exphcit forms ()B1|) and ()B2|) . The symbol A denotes the lattice 
Laplacian (with unit lattice spacing): 

The corresponding population-level description can be obtained from ()B6|) by letting N oo 
which eliminates the term of order A^~^ and allows us to replace {nj) by (nj)^, as in section 
im This leads to an equation for (pi = {rii) /N which is given by 

^ = mA0i - + 2h(j)i {I -(pi)- d(Pi . (B8) 
dT 

The final step that has to be taken in order to make contact with the equations used in the 
traditional approach, is to move from the lattice to the continuum. To do this we need to 
introduce a lattice spacing of e and take it to zero so that 

where the lattice site i is now replaced by the position vector x. In addition, the migration 
parameter has to be redefined, in order to absorb a factor of e^. The resulting equation is 
fl37|) . given in the main text. 

The derivation of the population-level description for the second version of the one- 
species model follows similar lines. The particular differences between this version and the 
one discussed above are described in the main text, and specifically by the equations ()39|) 
and pn|). The transition probabilities for this second version are 

T[. . .rii-l.rij . . .\. . . rii, rij . . .) = — riin 



T{. . .rii + l.rij . . .\. . .rii.rij . . .) = ^ (1 - 
T{. . .rii + l,nj - 1 . . . \ . . .ni,nj . . .) = (1 - ni)nj , (BIO) 

with similar equations with i and j interchanged, and 

T(. . . - 1 . . . I . . . . . .) = . (Bll) 

Note that the transition probabilities in ()B10|) are zero unless rij and rij are both equal to 
1 (competition) or rij = and rij = 1 (birth and migration), as required. The factors zN 
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and account for the choices of sites i and j and replace zQ and Q respectively in the first 
version. 

The master equation resembles ()B3|) . except that the single site processes are now re- 
stricted to the death process, and the two-site processes are more extensive, involving birth, 
competition, and migration. Defining 

we find using ()B4j) and the decoupling approximation {niUj) = {ni){nj), 

^{m) =mA{n,) - c{n,) { - > ]{n,) |> + 26 (1 - {n,)) { - > ]{nj) } - d{n,) . (B13) 



dr 



Denoting (jii) as 0j, the terms in curly brackets become 0(x, t) in the continuum limit, and 
so once again we recover equation ()37|) . given in the main text. 

The description of the IBMs when two species are present parallels that for one species. 
This similarity also holds for the initial stages of the derivation of the population-based 
equations, and so our description will be brief for both of these aspects. 

For the first version of the model, the transition probabilities for birth, competition and 
death processes are generalizations of ()13|) (the modifications are exactly the same as those 
made on © to give ()B1|) ). Those for migration of A's are ()B2|) . but with m replaced by nii 
and — Uj replaced by — rij — rrij (or N — rii replaced by A^ — rij — rrii). For migration 
of i?'s, they have the same form, but with mi replaced by m2 and with the substitutions 
rii rrii and nj ^ rrij. The master equation for P{n,m,t) is as before, but now including 
the greater number of allowed processes. There are two rate equations, found by substituting 
the master equations into 

d{nk) v^v^ dP{n,rh,t) d{mk) dP{n,m,t) 

{n} {m} {n} {m} 

Defining new quantities 

c., = ^^, &« = ]^, da = , ^a = ^, r = -, (B15) 



as in ()B5|) . we now let A^ ^ oo and replace averages of products by products of averages to 
obtain the equations 

^ = rhiA(j)i + ^ {(pitpj - (pjilji) - cii(t)1 - 2ci20j^/'i + 26i0i {I- (pi- ipi) - dicpi (B16) 
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and 

dip. 



2mo 



1712 Alpi^ 

ar z 



^ {i}i<\)j - ijj(f)i)-C22i'i-2c2iA(l)i + 2b2A {l-cpi- tjji)-d2^i . (B17) 



Here (pi = {ni)/N and tpi = {mi)/N. Writing (pitpj - (pjtpi as - tpi) - ijji{(f)j - (pi) we 

obtain equations ()44p and ()45|) in the continuum limit. 

For the second version of the two species model, the transition probabilities are general- 
izations of the one-species forms given by ()B10|) and ()B11|) . Specifically, for the competition 
process the term cniUj becomes Cunifij and, in addition, there are transition probabilities 
which are proportional to Cunirrij , C2iminj and C22'^i^j- For the birth process, the factor 
(1 — Ui) is replaced by (1 — rij — mi), and huj by hiUj or h2mj. The same holds for migration, 
but with 6, hi and 62 replaced by m, mi and m2 respectively. Finally, for the death process 
dui is replaced by diUi or d2mi. The master equation is straightforward, but tedious, to 
write down. 

Defining new quantities 



Ca/3 — — 7^, On 



lj,ba ; _ (1 - 

(Xrv 



a ^ fim^ 



t 



AT ' N ' N ' iV ' ■ iV ' ^^^^^ 

we find using the decoupling approximation — in which averages of products of any two of 

the variables {n, fh} are replaced by the products of their averages — that 



d 

d^ 



{m) = miA{ni) 



2ci2(n,) <j ^ j + 2bi (1 - im) - {mi)) |^ " ' (^^9) 



and 

^{mi) = m2A{mi) + 
dr 



- C22 {m 



2c2i(mi) I - ^{rij) > + 262 (1 - {ui) - (mi)) < - ^{nij) I - . 



(B20) 



Defining (rij) and {m^) as (pi and tpi respectively, we recover equations (j44p and (j45p in the 
continuum limit, up to slightly different definitions of the birth, competition, migration, and 
death rates. 
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Table 1 : Parameter values for Figures 1 and 2. 
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Table 2 : Parameter values for Figures 3-5. 
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Table 3 : Parameter values for Figures 6-9. 
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Figure Captions 



Figure 1 A comparison of tlieory to simulation for a single patch containing a single species. 
The upper panel shows the time evolution of the population density = {n)/N, where 
n is the number of individuals and is the size of the patch. The lower panel shows the 
time evolution of the variance, v{t) = ((n^) — (n)^) /N^, of the population. The subscript A 
refers to the fact that the individuals belong to species A. In this case the patch size has the 
relatively large value of 100 and we see that theory and simulation are in good agreement. 
See Table 1 for specific parameter values used in Figures 1 and 2. 

Figure 2 The same as Figure 1 but with the patch size reduced to 10. In this case the mean 
population density falls to zero due to fluctuation induced extinctions. The true variance 
first exceeds the large prediction and then falls steeply below due to the growing number 
of realizations which have become extinct. 

Figure 3 Comparison of theory to simulation for a patch containing two species, A and B. 
Initially each species has a density of one quarter of the total patch capacity. The upper 
panel shows the time evolution of the densities of A and B {(f){t) and ip{t), respectively), 
while the lower panel shows the time evolution of the variance of A compared with the large- 
A^ theory. In this case the A and B individuals have identical birth, death, and competition 
parameters, and the interspecific competition is set to zero. The patch size has the large 
value of 400. See Table 2 for specific parameter values used in Figures 3-5. 

Figure 4 The same as Figure 3 but with a patch size of 200. Note that although the mean 
field theory still works fairly well, the large- A^ prediction for the variance is very poor. Even 
for such a large patch, the fluctuations are increasing steadily with time. 

Figure 5 A similar scenario to Figure 3, with the addition of weak interspecific interactions 
(one fifth of the strength of the intraspecific interactions, with A out-competing B). Here, 
the patch size is 400, and the densities are fairly well approximated by mean field theory, 
although we see a slow decline in the B population. However, the fluctuations are again 
increasing with time and are not well-described by the large- A^ theory for large times. 

Figure 6 Comparison of mean field theory (smooth lines) and simulation (erratic lines) for 
two species A and i? in a spatial setting in which initially the A individuals occupy the left 
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half of the system and the B individuals the right half. The upper (lower) panel shows early 
(late) times. Here, the patch size has the relatively large value of 100. A out-competes B 
(meaning that C21 > Cu) but has a higher death rate, and so is invaded by B. See Table 3 
for specific parameter values used in Figures 6-9. 

Figure 7 The same as Figure 6 but here the patch size is reduced from 100 to the relatively 
small value of 10. Note that mean field theory is still in fairly good agreement with the 
high density B population, but shows significant deviation for the stochastically weakened 
A population. 

Figure 8 A similar scenario to Figure 6, but now A and B have identical death rates and 
yet B has less intraspecific competition than A. In this example A invades B. The patch 
size here has the relatively large value of 50. 

Figure 9 The same as Figure 8 except that the patch size is reduced from 50 to 10. Note 
that the invasion of A into B is severely slowed due to the stochastic weakening of A. 
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